6.5 SDM联合建模

6.5.1 sdm_PACKAGE

## sdm存在缺点，无法调参：还是优先使用ssdm和biomod2；
## 剩余

## 构建多模型建模：
## 使用sdm：
## 学习使用sdm快速构建多物种模型：



### 常规加载数据形式：####
file <- system.file("external/pa_df.csv", package="sdm")
df <- read.csv(file)
head(df)
d <- sdmData(sp~b15+NDVI,train=df)
d
# or simply:
d <- sdmData(sp~.,train=df)
d
## 加载分类变量：####
file <- system.file("external/pa_df_with_xy.csv", package="sdm")
df <- read.csv(file)
head(df)
d <- sdmData(sp~b15+NDVI+categoric1+categoric2+coords(x+y),train=df)
d
# categoric1 and categoric2 are categorical variables (factors), if not sure the data.frame has
# them as factor, it can be specified in the formula:
d <- sdmData(sp~b15+NDVI+f(categoric1)+f(categoric2)+coords(x+y),train=df)
d
# more simple forms of the formula:
d <- sdmData(sp~.+coords(x+y),train=df)


## 构建多物种同时建模：####
file <- system.file("external/multi_pa_df.csv", package="sdm")
df <- read.csv(file)
head(df)
# in the following formula, spatial coordinates columns are specified, and the rest is asked to
# be detected by the function:
d <- sdmData(~.+coords(x+y),train=df)
d
#--- or it can be customized wich species and which covariates are needed:
d <- sdmData(sp1+sp2+sp3~b15+NDVI+f(categoric1) + coords(x+y),train=df)


## 构建预测变量集：####
file <- system.file("external/predictors.grd", package="sdm") # path to a raster object
r <- brick(file)
r # a RasterBrick object including 2 rasters (covariates)
plot(r)
# now, we can use the species points and predictor rasters in sdmData function:
d <- sdmData(sp1+sp2+sp3~b15+NDVI,train=p,predictors = r)


## 与qg为例：####
library(raster)
library(tidyverse )
library(dismo)
library(SDMtune)
library(sdm)

setwd("E:\\闲鱼项目\\14赤皮青冈\\")
# install.packages("tidyverse")
## 加载物种分布数据：
cpqg <- read.csv("E:\\闲鱼项目\\14赤皮青冈\\建模最终使用的物种分布数据.csv")
cpqg <- cpqg[,1:3]

## 方法2：E:\闲鱼项目\14赤皮青冈\建模数据\气候建模数据集
biox <- paste("E:/闲鱼项目/14赤皮青冈/建模数据/气候建模数据集/","wc2.0_bio_30s_",c("02","07","11","12","14","15"),".tif.asc",sep = "")
bios <- stack(biox)

## 按照sdm的套路：需要提取点值：
re_cpqg <- raster::extract(bios,cpqg[,2:3])
re2 <- cbind(cpqg[,2:3],re_cpqg) %>% na.omit(.)

bg_cpqg <- read.csv("./建模数据/建模最终使用的背景分布数据.csv")[,2:3]
bg_sd <- raster::extract(bios,bg_cpqg) 
bg2 <- cbind(bg_cpqg,bg_sd) %>% na.omit(.)

c1 <- rep(1,dim(re2)[[1]])
c2 <- rep(0,dim(bg2)[[1]])
sp <- c(c1,c2) %>% data.frame(.)
names(re2)[1:2] <- c("x","y")
names(sp) <- "sp"

df <- rbind(re2,bg2) %>% cbind(.,sp) %>% data.frame(.)
head(df)
names(df)

## 建模气候数据集asc形式：
biox <- paste("",
              "wc2.0_bio_30s_",c("02","07","11","12","14","15"),".tif.asc",sep = "")
bios <- raster::stack(biox)

setwd("E:/闲鱼项目/14赤皮青冈/建模数据/气候建模数据集/")
bio1 <- raster("./wc2.0_bio_30s_02.tif.asc")
bio2 <- raster("./wc2.0_bio_30s_07.tif.asc")
bioxx <- brick(bio1,bio2)

d <- sdmData(sp~.+coords(x+y),train=df)

m <- sdm(sp~.,data=d,methods=c('rf'),replicatin='boot',n=4)





# ensemble using weighted averaging based on AUC statistic:
p1 <- ensemble(m, newdata=bios, filename='ens.img',setting=list(method='weighted',stat='AUC'))
plot(p1)
# ensemble using weighted averaging based on TSS statistic
# and optimum threshold critesion 2 (i.e., Max(spe+sen)) :
p2 <- ensemble(m, newdata=preds, filename='ens2.img',setting=list(method='weighted',
                                                                  stat='TSS',opt=2))
plot(p2)

6.5.2 SSDM

## 针对每个区域，每个分布地物种构建二元分布图：
library(SSDM)
library(raster)
gui()
## 加载物种分布数据：
xh_na <- as.data.frame(read.csv("D:/XH/第二阶段/xh/na.csv")[,2:3])
xh_as <- as.data.frame(read.csv("D:/XH/第二阶段/xh/as3.csv")[,2:3])
xh_eu <- as.data.frame(read.csv("D:/XH/第二阶段/xh/eu.csv")[,2:3])
xh_au <- as.data.frame(read.csv("D:/XH/第二阶段/xh/au.csv")[,2:3])
xh_sa <- as.data.frame(read.csv("D:/XH/第二阶段/xh/sa.csv")[,2:3])
occs <- list(xh_sa,xh_as,xh_au,xh_na,xh_eu)
names(occs) <- c("sa","as","au","na","eu")

## 加载环境数据：asc；
## 南美洲：
sa_envs <- stack(list.files("D:/XH/fifth/sa_area/",pattern = "asc",full.names =T))
as_envs <- stack(list.files("D:/XH/fifth/as_area/",pattern = "asc",full.names =T))
eu_envs <- stack(list.files("D:/XH/fifth/eu_area/",pattern = "asc",full.names =T))
au_envs <- stack(list.files("D:/XH/fifth/au_area/",pattern = "asc",full.names =T))
na_envs <- stack(list.files("D:/XH/fifth/na_area/",pattern = "asc",full.names =T))

## 分布点可视化：  #####
## 加载高程地图：H盘下的那个可视化分布范围：

## 构架可视化掩膜范围：##先放在这里！
par(mfrow = c(2, 2), mar = c(2, 2, 3, 1))
sa_e <- extent(sa_envs$BIO12)
plot(sa_envs$BIO12, xlim=c(sa_e@xmin-1,sa_e@xmax-5), ylim=c(sa_e@ymin-1,sa_e@ymax+1),legend = FALSE)
points(xh_sa,pch =21,col="red",cex =1)

as_e <- extent(as_envs$BIO12)
plot(as_envs$BIO12, xlim=c(as_e@xmin-1,as_e@xmax-5), ylim=c(as_e@ymin-1,sa_e@ymax+1),legend = FALSE)
points(xh_sa,pch =21,col="red",cex =1)



## 构建enm：#####
## 两种想法：1是用调参后的sdm，另外一种是使用ensemble：
## 所有参数模型：
# 'all' allows to call directly all available algorithms. 
# Currently, available algorithms include Generalized linear model (GLM),
# Generalized additive model (GAM),
# Multivariate adaptive regression splines (MARS), 
# Generalized boosted regressions model (GBM),
# Classification tree analysis (CTA), 
# Random forest (RF), 
# Maximum entropy (MAXENT),
# Artificial neural network (ANN), 
# and Support vector machines (SVM).




## 构建背景点：

## 生成不包含研究区域的随机背景点：
bg.xy <- dismo::randomPoints(sa_envs$BIO12, 2000,p= xh_sa)
bg_xy <- as.data.frame(bg.xy,colnames(c("long","lat")))
names(bg_xy) <- c("longitude","latitude")
# plot(bg_xy$x,bg_xy$y,cex=0.5)
# points(xh_sa$longitude,xh_sa$latitude,col="red",cex =1,pch=21)

## 构建swd形式，提供存在-背景点的0-1数据；
pa <-c(rep(1,nrow(occs$sa)),rep(0,nrow(bg_xy)))
occ_sa <- cbind(rbind(occs$sa,bg_xy),pa)
## 在没有背景点情况下的，设置取样 
PA= list('nb'=2000,'strat'="disk")

## cv：提供数据交叉验证的方法：
## 其中PA表示物种在仅提供存着点时的取样策略；
## cv     character. Method of cross-validation used to evaluate the SDM (see details below).
## cv.param    numeric. Parameters associated to the method of cross-validation used to evaluate the SDM (see details below).
cv = 'k-fold'
## 前面是分割比例，后面是分割的重复次数；
cv.param = c(4,5)

## 评估二元图的矩阵：
## (True Skill Statistic) maximizes the sum of sensitivity and specificity
metric = 'TSS'
axes.metric = "prop.correct" ##参见perterson的论文；


## 选择auc指数高于0.75的变量参与最终建模；在单一模型中不需要使用；
select.metric = c("AUC") ## 可提供的选项还有-prop.correct，也可以同时算多个指数：
select.thresh = c(0.75)
## 各模型的参数设置：
## 以下默认即可：
# Generalized linear model (GLM)
# Generalized additive model (GAM)
# Multivariate adaptive regression splines (MARS)
# Generalized boosted regressions model (GBM)
# Classification tree analysis (CTA)
# Support vector machines (SVM)
# Maximum Entropy (MAXENT)
# Random Forest (RF)

head(occ_sa)
## 用于演示：构建建模所用所有模型的单独模型：

sa_enm_glm <- modelling('GLM',occs$sa, Env =sa_envs,Xcol = 'longitude',
                        Ycol = 'latitude',PA= list('nb'=2000,'strat'="random"),
                        cv = "k-fold",cv.param = c(4,5),metric = "TSS",axes.metric = "prop.correct" ,
                        verbose = FALSE)
sa_enm_maxent <- modelling('MAXENT',occs$sa, Env =sa_envs,Xcol = 'longitude',
                        Ycol = 'latitude',PA= list('nb'=2000,'strat'="random"),
                        cv = "k-fold",cv.param = c(4,5),metric = "TSS",axes.metric = "prop.correct" ,
                        verbose = FALSE)
sa_enm_rf <- modelling('RF',occs$sa, Env =sa_envs,Xcol = 'longitude',
                        Ycol = 'latitude',PA= list('nb'=2000,'strat'="random"),
                        cv = "k-fold",cv.param = c(4,5),metric = "TSS",axes.metric = "prop.correct" ,
                        verbose = FALSE)
sa_enm_svm <- modelling('SVM',occs$sa, Env =sa_envs,Xcol = 'longitude',
                        Ycol = 'latitude',PA= list('nb'=2000,'strat'="random"),
                        cv = "k-fold",cv.param = c(4,5),metric = "TSS",axes.metric = "prop.correct" ,
                        verbose = FALSE)

sa_enm_mars <- modelling('MARS',occs$sa, Env =sa_envs,Xcol = 'longitude',
                        Ycol = 'latitude',PA= list('nb'=2000,'strat'="random"),
                        cv = "k-fold",cv.param = c(4,5),metric = "TSS",axes.metric = "prop.correct" ,
                        verbose = FALSE)

tt <- sa_enm_mars@binary+ sa_enm_svm@binary + sa_enm_rf@binary + sa_enm_maxent@binary + sa_enm_glm @binary
plot(tt )
sa_enm_mars@evaluation
points(xh_sa$longitude,xh_sa$latitude)

## ENSEmble：


# dir.create("./sa_enm")
sa_ESDM <- ensemble_modelling(c('GLM', 'MARS','MAXENT','RF','SVM'), occs$sa,Env =sa_envs,
                            rep = 2, Xcol = 'longitude', Ycol = 'latitude',save = TRUE,path="./sa_enm",
                           cv = "k-fold",cv.param = c(4,5),PA= list('nb'=2000,'strat'="random"),
                           metric = "TSS",axes.metric = "prop.correct" ,
                           uncertainty = TRUE, tmp = FALSE,
                           ensemble.metric = c("AUC"), ensemble.thresh = c(0.75),
                           weight = TRUE, verbose = FALSE)


plot(sa_ESDM @uncertainty)
sa_ESDM@algorithm.evaluation
plot(sa_ESDM@binary)



sa_enm_glm_as <- modelling('GLM',occs$as, Env =as_envs,Xcol = 'longitude',
                        Ycol = 'latitude',PA= list('nb'=2000,'strat'="random"),
                        cv = "k-fold",cv.param = c(4,5),metric = "TSS",axes.metric = "prop.correct" ,
                        verbose = FALSE)
warnings() 
sa_enm_glm_as@evaluation



ESDM <- ensemble_modelling(c('CTA', 'MARS'), subset(Occurrences, Occurrences$SPECIES == 'elliptica'),
                           Env, rep = 1, Xcol = 'LONGITUDE', Ycol = 'LATITUDE',
                           ensemble.thresh = 0, verbose = FALSE)


plot(ESDM@projection, main = 'ESDM\nfor Cryptocarya elliptica\nwith CTA and MARS algorithms')
plot(ESDM@binary, main = 'ESDM\nfor Cryptocarya elliptica\nwith CTA and MARS algorithms')

## 
knitr::kable(ESDM@evaluation)

plot(SDM@projection, main = 'ESDM\nfor Cryptocarya elliptica\nwith CTA and MARS algorithms')
SDM_projection <- project(SDM,Env_new)